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We consider the finite temperature metal-insulator transition in the half filled paramagnetic Hub- 
bard model on the infinite dimensional Bethe lattice. A new method for calculating the Dynamical 
Mean Field Theory fixpoint surface in the phase diagram is presented and shown to be free from the 
convergence problems of standard forward recursion. The fixpoint equation is then analyzed using 
dynamical systems methods. On the fixpoint surface the eigenspectra of its Jacobian is used to 
characterize the hysteresis boundaries of the first order transition line and its second order critical 
end point. The critical point is shown to be a cusp catastrophe in the parameter space, opening a 
pitchfork bifurcation along the first order transition line, while the hysteresis boundaries are shown 
to be saddle-node bifurcations of two merging fixpoints. Using Landau theory the properties of the 
critical end point is determined and related to the critical eigenmode of the Jacobian. Our findings 
provide new insights into basic properties of this intensively studied transition. 

PACS numbers: 71.30.+h, 71.10.Fd, 71.27.+a 



I. INTRODUCTION 

The correlation driven Metal-Insulator Transition 
(MIT) at finite temperature, also known as the Mott 
transition, is today one of the most intensively studied 
phase transitions in solid state physics. The problem 
contains competing energy scales making it inaccessible 
to perturbative methods. The seminal work of Metzner 
and VollhardlP spurred a rapid development in this field 
by introducing the limit of infinite connectivity. In this 
limit the Dynamical Mean Field Theory (DMFT^H be- 
comes exact and the lattice problem can be mapped to an 
auxiliary impurity model connected to a non-interacting 
bath. 

The paramagnetic MIT of the Hubbard model on the 
Bethe lattice where spatial correlations and magnetic or- 
der parameters are suppressed, has been studied by many 
authors in this particular limit P^l Regarding this transi- 
tion driven by Hubbard repulsion U, the emerging con- 
sensus is that it is a first-order phase transition terminat- 
ing at a critical point. The low temperature first-order 
transition line U C (T) is surrounded by a hysteresis re- 
gion with borders U c i(T) and U C 2(T), containing both a 
metallic and an insulating solution. At a critical temper- 
ature T c , the lines U cl (T), U c2 (T) and U C (T) all meet in 
the second-order critical end point, (U, T) — (U C (T C ), T c ), 
as schematically shown in Fig. [T] 

A theoretical framework for the understanding of the 
critical point has been presented by Kotliar et. aZ.j^ 
explaining it in terms of a DMFT Landau functional and 
an emerging zeroth mode in the "fluctuation matrix" of 
this functional. Moreover, the double occupancy D act as 
the thermodynamic conjugate variable to the (external) 
field t/P 

In this paper we show, using the impurity solvers Exact 
Diagonalization (ED) and Iterated Perturbation Theory 
(IPT), that this theory also can explain the existence 
of a third thermodynamically unstable solution in the 
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FIG. 1. (Color online) Sketch of the (U,T) phase diagram 
showing the first order thermodynamic transition line U C (T) 
and the instability edges of the insulating and metallic solu- 
tions, Ud{T) and U C 2{T) respectively. 



hysteresis region, previously reported by Tong and co- 
workers.^ We also present a general algorithm for finding 
fixpoints to the DMFT equations that is free from the sta- 
bility and con verge nce problems encumbering both for- 
ward recursio n 1 7 ! 11 ! and Newton methods in the vicinity 
of the hysteresis boundaries and the critical point. The 
algorithm is quite general and can be implemented with 
any DMFT impurity solver. Furthermore we present a 
method to calculate the Jacobian of the DMFT recursion 
relation at a fixpoint in the framework of ED and IPT. 
The properties of the Jacobian is then used to explain 
the origin of the numerical problems of forward recur- 
sion and Newton methods and how these difficulties are 
avoided by our algorithm. 

This paper is organized as follows: In Section [TTJ we 
give an introduction to the single band Hubbard model 
on the Bethe lattice, in Section II A we introduce DMFT 
and how it can be reformulated as a fixpoint problem. In 
Section [II B| we present our implementation of the Exact 
Diagonalization impurity solver and how the Jacobian of 



2 



the DMFT fixpoint function is calculated in this context. 



Section II C is used to explain the same details for the 
Iterated Perturbation Theory impurity solver. Based on 
the general fixpoint problem we investigate the local con- 
vergence properties of the mentioned fixpoint solvers in 
Section |II D| In Section |II E| we give a brief description 
of the thermodynamics of the MIT. In Section |III| the 
results are presented and put into relation with previous 



work in Section IV Finally we give a short conclusion in 
Section [V] 



II. THEORY 

The Hamiltonian H of the half-filled Hubbard model 
is given by, 
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FIG. 2. (Color online) The self consistent DMFT equations 
and the two possible fixpoint function formulations Fup(El) 
and F uf j(G ). 



with nearest neighbor hopping —t, local Hubbard repul- 
sion U and the chemical potential, fx — U/2. In the limit 
of infinite dimensions, d — > oo, the hopping matrix el- 
ement t has to be rescaled as, t — > t/yd, in order to 
retain a finite kinetic energyP On the Bethe lattice the 
non-interacting density of states p(°) (uS) is semicircular 
and given by, 
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where W is the bandwidth (W = At), In this study we 
use, W = 2eV. Since the Bethe lattice is bipartite the 
ground state of H is antifcrromagnetic at low tempera- 
ture, which should in principle suppress the MIT of the 
paramagnetic state studied here. In the spirit of previ- 
ous studies^ we enforce the paramagnetic state by impos- 
ing translational invariance ignoring spatial correlations. 
This also enables us to connect to experimental results 
on more frustrated multi-band systems not displaying the 
antiferromagnetic instability. 



acting as a dynamic Weiss field. The bath Green's func- 
tion Go is obtained by subtracting the local interactions 
using E L . 

G (iuj n ) = [G7 1 («w„) + E L (zw„)] _1 (4) 

The local interactions of the lattice Hamiltonian H and 
Go now fully determine the action Si for the impurity 
system, 

rP 

Si[G ] = U / dr c\(T)c t (T)cl(T) H (T) - 
Jo 

- f dr f dr'£ct(r)Go V-r'Mr'). (5) 
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Solving the impurity problem is still a formidable task 
but the single impurity system is now well within reach 
for state of the art numerical algorithms. Using an impu- 
rity solver the paramagnetic interacting impurity Green's 
function, Gi(ioj n ), can be calculated as, 



A. Dynamical Mean Field Theory 

Let us first introduce some concepts used in the formu- 
lation of DMFT. In the limit of infinite coordination the 
lattice self-energy E^ is local, E^(R, iui n ) — > Ei(iw n ), 
and DMFT^EI [ s an exact theory. In the continuum limit 
the local lattice Green's function, Gl(R = 0,iuj n ) = 
G L (iu) n ), is given by, 



G L (iu n ) = duj - 
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where iu> n are Matsubara frequencies. With the lattice 
we associate an auxiliary impurity connected to a bath 



Gi a<T '(iLO n ) 



dre^{T ClT (T)cl,(0)} Sl[GQ] , 



where, G/ = G m = G IU , G m = G m = . (6) 

The corresponding impurity self-energy Ej is calculated 
by inverting the Dyson equation, 



E/(iw„) = G 1 (iu) n ) - Gj 1 {iuj n ) . 



(7) 



Now given E^, the equations, ([3j|4j[6]and[7|, can be used 
to calculate a corresponding E/. Incorporating these 
steps into a single DMFT function F gives the short form, 



F ufi (T. L ) = Y. I , 
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FIG. 3. (Color online) Truncated SIAM, with one correlated 
impurity level ei, bath levels e*., 1 < k < Nb, hybridizations 
Vk and Hubbard repulsion U. 



where the Hubbard U and the inverse temperature, /3 = 
1/ksT, are external parameters. 

The last step is to find self consistent solutions where 
the lattice and impurity self-energies coincide, Ei(iu; n ) = 
Tii{iui n ). This is equivalent to finding fixpoint solutions 
E* of the DMFT function F Uj3 , 



(9) 



In principle any fixpoint algorithm for multidimensional 
functions can now be applied on F\jp to find DMFT solu- 
tions. The choice of E as the fundamental variable in the 
fixpoint scheme is not unique, as one can alternatively use 
the bath Green's function Go as fixpoint variable giving 
fixpoint solutions, Fjjp(Gq) — Gq. Both possibilities are 
indicated in Fig. [2j where the coupled equations forming 
the DMFT fixpoint functions are shown schematically. 



B. Exact Diagonalization 

To solve the impurity problem of Eqs. ^ and (|6| we 
have implemented^ the Exact Diagonalization (ED) al- 
gorithm by Caffarel and Krauth.^ In ED the impurity 
problem is projected to a truncated Single Impurity An- 
derson Model (SIAM), with the Hamiltonian, 

H S1AM = }Xei - p)c\c a + ^ e kcl a c k a+ 

o ka 
ka 

keeping the local interaction of the lattice Hamiltonian 
H but replacing the hopping term with hybridization Vk 
to a set of "bath" states at energies e k , see Fig. [3j 

As H is particle-hole symmetric the same symmetry is 
imposed on the SIAM parameters. Letting e/ = 0, the 
bath level energies efc are placed symmetrically around 
zero energy. Moreover for each pair, ek = — e^. 7^ 0, 
the hybridizations are the same, Vk = Vz. Half filling is 
obtained by fixing the chemical potential [i to, fi = U/2. 



The first step of the algorithm is to project the im- 
purity bath Green's function Go to the non-interacting 
(U = 0) SIAM Green's function Gg IAM , 
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through minimization of a penalty function x 2 , 
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with respect to the SIAM parameters e k and To 
minimize \ 2 a standard conjugate gradient minimization 
algorithm is used. Many forms of the penalty function 

x 2 



be constructed^! and in this work we choose Eq. 



(12 1, which is sensitive to the low frequency behavior of 
Gq iam , and N = 2 9 matsubara frequencies. 

With the parameters of H S1AM determined its matrix 
representation in the occupation number basis is calcu- 
lated. The symmetries of H S1AM can be used to block- 
diagonalize the matrix representation. For simplicity 
we only exploit the symmetry that subspaces with fixed 
number of spin up and spin down are not mixed 
by H s iam ■ Let us denote the eigenstates of Hsiam by \ y ) 
where, H slAM \v) = E v \v). 

The eigenstates are explicitly calculated by diagonal- 
ization and used to calculate the impurity Green's func- 
tion G/o-o-' from the Lehmann spectral representation,^ 
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(13) 

As the the system is paramagnetic, Eq. ([6| holds, and 
only one impurity Green's function Gi has to be calcu- 
lated. The double occupancy D is given by, 
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The described steps of the ED algorithm are schemat- 



ically shown in Fig. |4jand replaces the, Gj = (cc^)si[G a ]> 
block in the DMFT equations of Fig. [2] 

In our calculations the dimension of the SIAM Hilbert 
space is the limiting factor of the ED algorithm. With 
one impurity level the number of fermionic states Nf is 
given by, Nf = Nf, + 1, where Nb is the number of bath 
levels. The corresponding size of the Hilbert space be- 
comes 2 2N J , growing exponentially with respect to Nf. 
Dividing the Hamiltonian in blocks of constant (n^,n^) 
gives a set of smaller Hilbert spaces with dimensions, 
(n / )(rv / ) - Ovrr calculations converge rapidly with the 
number of fermionic levels Nf and for, Nf = 6, used in 
our calculatio ns, th e region around the critical point is 
well converged! 16 * 17 ^ 



4 



Go 







r 


- 


Gi = 


cc t ) S/ 






Gi 





[e k ,V k ] 



> f 



1 1/) 



G/ = E 



FIG. 4. (Color online) Detailed schematic (right) of the Exact 
Diagonalization impurity solver (left). 



1. Fixpovnt function and J acobian 

We now discuss the coupled DMFT equations as a fix- 
point function in terms of Go instead of S^. In the ED 
algorithm Go is parametrized by a small number of pa- 
rameters, x = [efc,Vfc], when projected on to the SIAM 
and a solution of the DMFT equations can be formulated 
as a fixpoint problem in x, Fjjp(yL*) = x*. 

The reduction of parameters using x instead of 
Go(iui n ) facilitates a direct calculation of the Jacobian, 
Jf(x) = Vi r {/ I a(x). At each obtained fixpoint x*, 
Fupix*) — x* , the low dimensional parameter space al- 
lows us to use a modified central finite differences formula 
to calculate the Jacobian, 



JV(x) 



/t/^x + hk n ) - F V p(x - fox.,-, 
* 2h 



(15) 



where h is the discretization, x„ is the unit vector in the 
n:th dimension. Due to the large parameter spread in 
x* relative scaling was applied to stabilize the numerical 
evaluation of Jp (x*). 



C. Iterated Perturbation Theory 



£(tj) (or Green's function) using a finite number of Mat- 



subara frequencies, u) n 



2 f(n 



|), with, < n < 



N. Thus defining a "Matsubara-periodized" self-energy 
T,(iuj n ) through, E(rj) 

W \ TT^N-l i 

E(lCJ n ) = ^ 3= o e 



with period ^S-N. We formulate a fixpoint equation in 



/3 _1 £„=o , e-^E^n), and, 
Tj T,(Tj), which is now periodic 
La- 
terms of the finite dimensional self-energy similarly, 

Fup (E(iw„)) = E'(iwn) , (16) 

allowing for the same study of the fluctuations around 
the fixpoint as for the ED formulation but now in the 
TV dimensional space spanned by E(ia; n ). Specifically we 
calculate the Jacobian of Fup through, 



Fup{Y, + hz n ) - F Uf j(T, - hz n ) 
2h 



(17) 



where, z n = i(x n — X- n -i)/y/2, is a unit vector in the 
particle-hole symmetric £(iw„) subspace and h is the fi- 
nite difference discretization. 

We will not present the details of the calculations here 
in terms of "Matsubara-periodized" Greens functions^ 
but only point out some main features. Since we have 
discretized the Green's function, we have to be particu- 
lar careful about defining, G(r 3 = 0). Particle-hole sym- 
metry implies that, Gl(t) = — Gl(— r), and to preserve 
this we define, G(ij = 0) = 0. With this definition the 
first-order Hartree-contribution to E is zero and corre- 
spondingly, /j, = 0, at half-filling. The discretized IPT 
approximation for the self-energy is given by, E(Tj) = 
-[WtTjOGt-Tj) = U 2 G 3 (Tj), with G(Tj = 0) = by 
definition. From this follows also that E is purely imagi- 
nary. 

The DMFT equations can be written exactly in terms 
of the periodized Green's function and self-energy by re- 
placing Eq. ([3| with, 



G L {iuj n ) 
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P°M 

(coth (iu) n - w))" 1 - S(iw n ) 



(18) 

where the coth(...) term is the exact expression for the 
discrete Fourier transform of the non-interacting Green's 
function on the lattice. 



To show the generality of the fixpoint analysis we also 
consider the Iterated Perturbation Theory (IPT/p for- 
mulation of the DMFT equations, which amounts to 
solving the impurity problem, Eq. ([6]), perturbatively 
to second order in U. To discretize the problem in 
this case we consider directly the iteration scheme of 
the self-energy E(r) defined on a discrete time mesh, 
tj — with N a constant integer and j integer. Clearly 
the step size (3/N needs to be increased with decreas- 
ing temperature to be able to capture the self-energy 
or Green's function sufficiently well. Because of the 
discretization in time, we can represent the self-energy 



D. Fixpoint solvers 

A common ingredient in all DMFT calculations is 
solving a fixpoint problem. Conceptually the choice to 
parametrize the fixpoint function in terms of E or Go 
is irrelevant, since the resulting fixpoints are equivalent 
whichever coordinates are used. In this section we adapt 
E as the fixpoint-variable keeping in mind that it is ex- 
changeable with Go (or x = [e^., V&]). 

We now discuss the two most widely used algorithms 
for solving the fixpoint problem, forward recursion and 
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Newton methods and for each method the local conver- 
gence properties around a fixpoint £* will be explained 
in terms of the dominating eigenvalue e of the Jacobian, 
Ji?(S*). Finally we introduce the phase space extension 
and explain why this method is free from some of the 
deficiencies of forward recursion and Newton methods. 



1. Forward recursion 

The most common algorithm for solving the DMFT 
equations is the fixpoint forward recursion^ Given an ini- 
tial guess S a series {En} is generated by the recursion 
relation, 

Fupp n ), (19) 

and a fixpoint £* is found if the series converges, E* = 
Eqo. To study the convergence properties of the series 
{£„} in the vicinity of a fixpoint £*, we can approximate 
Fjjp by its first order Taylor expansion, 

F ufi (£* + ST) « F up (£*) + J F (£*) ■ ST , (20) 

where Jf is the Jacobian matrix of Fup, Jf(T) = 
Vi r [/ i a(E) and ST is a small perturbation. By repeated 
application of the recursion near the fixpoint, 

E =£* + ST 

E„ =£* + J F (E*) n <SE, (21) 

we easily observe that the convergence is determined by 
the eigenvalues of Jf and in particular by the eigenvalue, 
e, with the largest magnitude. Hence we require, |e| < 1, 
for forward recursion to converge at all. This imposes 
a restriction on the solution space that can be found by 
this scheme. 

If the Jacobian at a fixpoint has an eigenvalue larger 
than one in magnitude, the algorithm will only converge 
if the perturbation ST has no components in the corre- 
sponding eigenspace. Any contribution in ST from this 
eigenspace will be amplified and T n will move away from 
the fixpoint £*, and the forward recursion algorithm will 
be unable to find the fixpoint in the first place. Fixpoint 
forward recursion can therefore only be used to find a 
subset of all fixpoints T* of the function Fup whose Ja- 
cobian have all eigenvalues bounded by one in magnitude. 
And if, |e| — » 1 — , when tuning an external parameter the 
fixpoint forward recursion experiences a critical slowing 
down of convergence, due to the damping factor e™ of a 
perturbation. This phenomena has been reported for the 
MIT of DMFT when appr oach ing the hysteresis bound- 
aries of the phase diagram!^^ 

2. Newton methods 

The family of Newton's method and the quasi New- 
ton methods^ are all multi dimensional root solvers with 



better stability properties than forward recursion. Broy- 
den's method^ from this class of algorithms have re- 
cently been applied in the context of DMFT.I^il In order 
to use a root solver, the fixpoint problem in Eq. (|9| is 
simply reformulated to a root problem, 

R up (T*)=F ul3 {T*)-T* =0, (22) 

where the Jacobian Jr of Rup has the form, 

J R (T) = VR V p{T) = VF uf) {T) - 1 = J F (T) - 1 . (23) 

The series {£«} is in the case of Newton's method 
generated as, 

£„+! = T n - {MVn^Rup&n) ■ (24) 

In the linear regime in the vicinity of a fixpoint E* 
where, R V p{T* + ST) w J R {T*)8T and J R {T* + ST) « 
Jr(E*), the series converges in one iteration as, 

E =T* +5T 

Ex =E - J fl 1 (So)i?c//3(E )«E*, (25) 

assuming that Jr(E*) is invertible. Translated to the 
eigenvalue spectrum of the Jacobian Jp{T*) of the fix- 
point function Fup(E*), local convergence is achieved as 
long as no eigenvalue is equal to one, a less restrictive 
requirement compared to forward recursion. Newton's 
method will therefore converge even in areas of parame- 
ter space where forward recursion fails completely. 

3. The phase space extension 

Although Newton's metod allows us to work with, 
|e| > 1, a problem remains when, |e| = 1, which (as we 
will show) is precisely at the hysteresis boundaries of the 
MIT. To be able to trace the solutions across this singu- 
larity we reformulate the problem in an extended phase 
space where the resulting Jacobian no longer becomes 
singular at the hysteresis boundaries. 

We construct a real- valued function A(T) whose value 
lifts the degeneracy of the coexisting fixpoints of the 
DMFT fixpoint function. Using A(T) we write down 
an extended root problem that not only finds a DMFT 
solution but also fixes the value of A to some given value 
a, A(T) = a. By increasing the dimension of the root 
problem by one, the extended root function R a p can be 
defined as, 

R a p(U,T) = (a-A(T),Rup(T))=0, (26) 

where U is treated as a free parameter to be varied along 
with E in order to obtain, a — A(T) — 0. By constraining 
the value of A(T) to a the degeneracy of the fixpoints is 
lifted and the resulting Jacobian is invertible even at 
the hysteresis boundaries. 

In this investigation, A(T) — Im[E(iwo)], was used and 
sampled on a logarithmic grid. The extended root prob- 
lem was solved using Broyden's second methodP^ As 
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Broyden's method does not require any Jacobian eval- 
uations it has the same computational cost as forward 
recursion, but with an almost constant convergence rate 
in the entire phase diagram. 



E. Thermodynamics 

Studying the MIT in the (U, T) phase space requires 
understanding of the thermodynamics of H in Eq. ([lj. 
Let us first consider the expectation value (H — fj,N) and 
introduce some notation, 



(H - fiN) =(f) +U^2(D i ) - A ^(n JCT ) 



= (T> + UN (D) 



(27) 



(28) 



where N is the number of sites, Di = n^n^ is the on 
site double occupancy, and the single particle hopping is 
contained in the kinetic term T. In the last step we have 
assumed half- filling with (n^) = 1/2 and /i = {7/2. It is 
now evident that U acts as an external field and can be 
considered a conjugate variable to the double occupancy, 
D = (D). From this it is possible to derive a Maxwell 
construction in U and D analogous to the formulation in 
density and pressure for the van der Waals equation of 
state.™ 

From the definition of the grand partition function 
Z{j3, (j,, U) and the free energy Q(/3, fj,, U), 



Z = e 



= Tr 



the derivative of £1 with respect to U is given by, 



on 

dU 







ld_ 



bxZ = N ( (D) - 



1 



(29) 



(30) 



The free energy difference AO between two points on an 
isotherm can be expressed as, 

dU [(D)-~) . (31) 



provided that there is an adiabatic connection between 
U\ and C/ 2 - 

In the case of the MIT this can be used to determine 
the thermodynamic first order transition given by the 
three DMFT solutions on an isotherm in the hystere- 
sis region, as reported previously by Tong el. This 
third unstable solution connects the metallic and insu- 
lating solutions making it possible to calculate the free 
energy difference between the metal and insulator for a 
given Hubbard U. 



III. RESULTS 

Using the phase space extension and Broyden's second 
method to solve the DMFT equations we find three so- 
lutions in the hysteresis region at fixed U and (3. One 
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FIG. 5. (Color online) Isotherm at, T = 22.50 meV, inside the 
hysteresis region, with the double occupancy D as a function 
of U . Shaded areas show the Maxwell construction. The solid 
line is composed of dense DMFT-ED solution points. 



metallic solution with high double occupancy D, one in- 
sulating solution with low D and a third intermediate 
"unstable" solution, see Fig. [5] As a function of U these 
three solutions form a continuous Z-shaped isotherm in 
D(U), where the unstable solution adiabatically connects 
the metallic and insulating solutions. The unstable solu- 
tion is not an artifact due to the phase space extension 
since it also is a solution to the DMFT root problem, 
Rup(E*) = 0, of equation Q. 



With the continuous isotherm D(U), of Fig. [5j it is 
possible to apply the Maxwell construction, Eq. ( pT| ), 
to determine the thermodynamic first order transition 
U C (T). Where the free energies of the metallic and insu- 
lating solutions coincide. This corresponds to equating 
the enclosed areas left and right of U c , as shown in Fig. 
[5] It is evident that the unstable solution always has a 
free energy higher than both the metallic and insulating 
solution, thus always being thermodynamically unstable. 
Increasing the temperature towards the critical temper- 
ature T — > T~ shrinks the size of the hysteresis region 
and at T = T c it disappear. 

Many of our results can now bee seen to be quite gen- 
eral properties of the MIT. The DMFT fixpoint solutions 
form a continuous surface in the (U, T, D) phase space. 
At the critical point {U c ,T Cl D c ) this surface has a cusp 
singularity. As a function of the DMFT recursion this 
surface opens up a pitchfork bifurcatiorP^l in the (T, D) 
plane and two of the solutions annihilate by saddle-node 
bifurcations^ at the hysteresis boundaries U C \(T) and 
U C 2{T) in the (U,T) plane. Studying the isotherms of 
D(U) and e(U), where e is the in magnitude largest eigen- 
value of Jp(x*), around the critical point, see Fig. |6j we 
can explain the behavior of the common algorithms used 
to solve the DMFT equations. 

Above the critical temperature, T > T c , the eigen- 
value e is always less than one and both forward recur- 
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FIG. 6. (Color online) Double occupancy D (upper panel) 
and maximum eigenvalue e of Jf(E*) (lower panel) plotted 
against U on isotherms above, close to and below the critical 
point (circles, triangles and squares respectively). Saddle- 
node bifurcation boundaries U c i and U C 2 occur when e — 1. 
The solid lines are composed of dense DMFT-ED solution 
points. 



sion and Newton's method converge. In this regime the 
maximum of e at the coupling U e deter mine s the center 
of the thermodynamic crossover regionP^H At the crit- 
ical temperature, T ~ T c , forward recursion displays a 
critical slowing down of convergence, as e — > 1~ when 
U — s- U c , while Newton's method becomes unstable first 
at the critical point (U C ,T C ). Below the critical tem- 
perature, T < T c , the behavior of e(U) becomes more 
complicated. Following the metallic solution (high D) 
the solution annihilate with the unstable solution at the 
second hysteresis boundary £/ c2 through a saddle-node 
bifurcation. This coincides with e — > 1~ and explains 
the critical slowing down of forward recursion when ap- 
proaching the hysteresis boundary from the inside of the 
hysteresis region. The behavior of the insulating solution 
(low D) at the first hysteresis boundary U c \ is analogous 
to that of the metallic solution. 

The unstable solution emerges at the hysteresis bound- 
aries U c i and U C 2 through the saddle-node bifurcations 
and has e > 1 in the entire hysteresis region. Thus for- 
ward recursion will never find this solution, although we 
have been able to trace it using Newton's method by 
supplying a close enough initial guess. But Newton's 
method is still unstable on the hysteresis boundaries and 
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FIG. 7. (Color online) Phase diagram in U and T plane, 
markers are ED-DMFT data and the dotted line correspond 
to, ao(u,t) = 0, to linear order i.e., /3ow + 7oi = 0. (Solid lines 
are guides for the eye.) 



a far better method is the phase space extension which 
converges everywhere in the studied parameter range. 

Recalling the definition of e(U, T) as the maximum 
eigenvalue of Jp(S*) let us dehne U £ (T) to be the cou- 
pling that maximizes e for fixed T, see Fig. [6] This allows 
for a precise determination of the critical temperature T c 
and critical coupling U C (T C ) as, e(U e (T c ), T c ) — 1 and 
U C (T C ) = U C {T C ). By linear interpolation of e{U e (T),T) 
from isotherms above and below T c the critical end point 
can be directly determined, as indicated in the lower 
panel of Fig. [6] Away from T c , U e {T) is not equal to the 
coupling U C (T) where the first order transition occur. 

The (U, T) phase diagram can now be represented by 
four lines, see Fig. [7| the hysteresis boundaries U C \{T) 
and U C 2(T) that confine the region of fixpoint coex- 
istence, the first order thermodynamic transition line 
U C (T) between metal and insulator and the maximal 
eigenmode curve U e (T) . Where U e (T) gives the crossover 
between metal and insulator for temperatures above the 
critical temperature T c . 

To accurately determine the location (U C (T C ),T C , D c ) 
of the critical end point and its critical properties we fit a 
Landau functional model L to the calculated DMFT-ED 
isotherms near the critical end point, 



d 2 d 3 
C(u, t, d) = a d + a\ — + a 2 — + 



2 

T c , d = 
a„(u,t) 



3 
D 



4 ' 

D c and a n are 



(32) 



where, u = U - U c , t = T 

linear functions in u and t, a n (u,t) — f3 n u + j n t. The 
expansion to fourth order in d is the minimal model for 
a system with a cusp singularity.^ The free parameters 
of the fit are, U c , T c , D c , /?„ and j n . 

For fixed u and t the extremal points of C, 



ddC — ciq + aid + a 2 d 2 



0, 



(33) 



are the Landau analogue of the DMFT fixpoints spanning 
a continuous surface Sc = {(u,t,d) : ddC(u,t,d) = 0}, 
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FIG. 8. (Color online) Landau fit (solid lines) and DMFT- 
ED data (markers) in the (u — r fo/Pot,T) plane close to the 
critical point. 



in (w, i, d) phase space. The hysteresis boundaries on Sc 
satisfy, d\L — 0, and in addition the critical point is 
characterized by, d\L = 0, d\L > 0. 

To obtain a fit to the DMFT-ED data, on the ±1 meV 
scale around the critical point, a second order tempera- 
ture term was added to ao(u,t), 



a (u, t) = j3 u + 7 i + j^t 2 • 



(34) 



Note that this extra term does not change critical behav- 
ior of the Landau functional, but changes the behavior 
of its "unstable" solution significantly. 

We asses the ability of the minimal cusp singularity 
Landau functional to describe the DMFT fixpoints in 
the vicinity of the critical point by comparing, the (U, T) 
phase diagram, the isotherms of D(U) and e(U) and fi- 
nally the critical behavior along the first order transition 
line. 

When studying the phase diagram in the (£/, T) plane 
it is useful to compensate for the linear slope of the 
hysteresis region at the critical point. From the Lan- 
dau functional the slope is obtained as the dotted line, 
(3qu + jot = 0, shown in Fig. [7] By the transformation, 

u — » u — jP-t, this line becomes vertical and the corre- 

r-i 

sponding transformed phase diagram is shown in Fig. [8] 
The phase diagram of the DMFT-ED data and the Lan- 
dau functional converge when approaching the critical 
point. 

The Landau functional isotherms of D(U) agree re- 
markably well with the calculated DMFT-ED data as 
shown in Fig. [9j An important additional fact is that 
the largest eigenmode e of the Jacobian Jp of the DMFT 
fixpoint function Fjjp and the second derivative of the 
Landau functional d 2 d Z are related through, 



djC = C(e - 1) 



(35) 



for, (u,t,d) £ Sc, where C is a constant factor. (See 
the lower panel of Fig. [9]) As only one eigenmode of the 
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FIG. 9. (Color online) Landau fit (solid lines) and DMFT- 
ED data (markers) in the vicinity of the critical point for the 
shifted double occupancy, d — D — D c (upper panel), and 
maximum eigenvalue e of JV(£*) (markers) and the scaled 
and shifted second order derivative of the Landau functional 
d^C/C+l (solid lines) (lower panel), plotted against, u = U — 
U c , on isotherms with, T = 25.0 meV + (0.650, 0.600, 0.575, 
0.550, 0.500) meV, for circles, triangles, squares, pentagons 
and hexagons respectively. 



DMFT Jacobian Jp(T,*) becomes critical in the hystere- 
sis region, this mode alone governs the critical behav- 
ior of D, while all other eigenmodes gives the universal 
structure of the phase diagram around the critical point. 
The observed proportionality of Eq. ( 35 ) confirms that 



our Landau functional, parametrized by the single pa- 
rameter d, is able to describe the final "effective" critical 
behavior. 



From Eq. (35) we can construct Ug2 C (T) analogously 
to U e (T) as the coupling maximizing d\L on an isotherm 
with temperature T. Comparing Ug2 C (T) to UJT), see 
Fig. |8j qualitative agreement is achieved below T c and 
quantitative agreement above T c . 

As a final test of the Landau fit we study the critical be- 
havior of the fixpoints on the surfaces defined by the first 
order transition line U C (T) and the critical eigenmode 
maximum line U e (T). On these surfaces, (U C (T),T, D) 
and (U e (T),T, D), the emergence of the pitchfork bifur- 
cation at T c is clear, and the Landau model and DMFT- 



ED data are again in good agreement. (See Fig. 10 ) 

We have shown that the Landau functional C quanti- 
tatively reproduces the critical properties of the DMFT 
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FIG. 10. (Color online) Double occupancy D as a function 
of temperature T on the {U C (T), T, D) and (U e (T),T, D) sur- 
faces from DMFT-ED (markers) and the Landau fit (solid 
lines). (Inset) Logarithmic plot of the metallic branch of 
(U C (T), T, D) showing the d ~ 1 2 critical exponent. 



fixpoint surface in the (U, T, D) phase space and criti- 
cal exponents can now be derived directly from £. For, 
U = U c , the critical behavior of the double occupancy 
with respect to T is given by Eq. (33 1 as, 



D -D r 



(36) 



and equivalently for T = T c the double occupancy as a 
function of U has the same critical form, 



D - D r 



\U-U e \*. 



(37) 



Along the first order transition line we have a very 
different critical behavior. Approaching criticality this 
line coincides with the first order term in the Landau 
functional being zero, ao(u, t) — /3qu + 7o£ + Jot 2 = 0. 
Sufficiently close to the critical point only the linear order 
contribute giving, u = ~^t- With this constraint on u 
we regain "Ising" scaling exponents and the temperature 
critical exponent changes to 1/2 ie., 



D - D r 



IT - T r 



iff U = U c 



To 
A) 



(T-T c ), (38) 



which is also confirmed by the DMFT-ED data and Lan- 
dau fit in the inset of Fig. 
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The Landau fit also gives a very precise location of the 
critical point, in the approximation of DMFT-ED with, 
N f = 6, we obtain D c = 0.03244 ± 0.0001 pairs/site, for 
U c and T c see Table [i] The error in T c is estimated by 
the temperature difference between the isotherms closest 
to T c and the error in U c and D c are estimated from 
the isotherms in Fig. [9] The values are compatible with 
previous reports on the critical point, see Table ||J Even 
though the precision is very high, the errors are estimated 
within the approximation of ED with a SIAM size of, 



TABLE I. Comparison of (U C ,T C , D c ) for the second order 
critical end point. 





u c \y v J 


T 

J- c 


ymc v I 


This work, ED 


2.3398 ± 0.003(PJ 


25.5625 


±0.0124^] 


HF-QMO 3 


2.3325 ±0.015 


27.5 


±0.2 


HF-QMCP 


2.38 ±0.02 


25.0 


±3.0 


eePI 


2.34 


25 




nrgP 




40 




This work, IPT 


2.46073±0.00050 


46.9048 


± 0.0550 




2.46315 


46.895 




ipqiizi 


2.51 


44.0 





Errors given with respect to ED using Nf 
finite size effects. 



: 6, not including 



12 The real 



he SIAM and 



Nf = 6, using the weight function of Eq. 
accuracy is lower due to the finite size of t 
the particular choice of weight function. 

The fixpoint surface Sc can be intuitively understood 
in terms of stationary points of the Landau functional C 
as a function of d. In Fig. [TTJ C is shown for fixed T 
and U on the hysteresis boundaries and at the first order 
transition. For U = U c the fixpoints correspond to two 
equal local minima and one unstable maxima. As the 
local minimas have the same value of C the system is un- 
stable and can undergo a first order phase transition. At 
the hysteresis boundaries U = U c i , U c % a saddle-node bi- 
furcation occur through the appearance of an inflection 
point, with two coinciding stationary points that sepa- 
rate when moving in to the hysteresis region. The local 
maxima corresponding to the unstable solution is always 
the stationary point with the highest free energy and can 
never be the thermodynamic ground-state. 

To demonstrate generality of the phase space extension 
and the critical properties of the MIT the same calcula- 
tions have also been performed using the impurity solver 
IPT. In the IPT calculations the Matsubara formalism is 
implemented in terms of periodizcd Green's functions. In 
the ED calculations the DMFT problem Fjjp was formu- 
lated as a fixpoint problem in terms of the parametrized 
bath Green's function Go(iuj n ), Fufs(Go) = Fy^(x), 
while in the IPT calculations the fixpoint problem was 
formulated in terms of the self-energy S(iw rl ), Fup(Y,). 
The Jacobian was evaluated numerically and the maxi- 
mal eigenvalue e determined. Applying the phase space 
extension to IPT we get exactly the same critical be- 
havior as before, see Fig. [l2j although the loc ation of 
the critical point is shifted as previously reported j^ZI see 
Tab. HI 



IV. DISCUSSION 

Using the phase space extension and the Jacobian 
Jf(£*) of the DMFT fixpoint function we have mapped 
out the metal insulator phase digram, the different 
DMFT fixpoints an d th eir bifurcations. This extends 
previous calculations^^ by describing the unstable solu- 
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FIG. 11. (Color online) The Landau function £ as a function 
of the double occupancy, d — D — D c , at T = 25.5 meV, 
slightly below T c , for U fixed at U&, U C 2 and U c - Markers 
indicate stationary points of C. 
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tion, which allows us to compute a continuous two dimen- 
sional fixpoint surface in the three dimensional (U, T, D) 
phase space. 

The critical point has been classified as a cusp catastro- 
phe with a pitchfork bifurcation on the first order tran- 
sition line. The hysteresis boundaries have been shown 
to be saddle-node bifurcations of two merging fixpoints, 
one stable and one unstable respectively. Using the ex- 



plicit calculation of Jp(S*) we have confirmed the pre- 
diction that a single critical eigenmode governs the MIT 
in DMFl^ and shown that this mode becomes critical 
not only at the critical point but also on the hystere- 
sis boundaries. The calculated fixpoint surface was then 
used to fit a cusp catastrophe minimal Landau functional 
expansion in the vicinity of the critical point. 

The excellent agreement between the mean field Lan- 
dau model and the DMFT-ED data shows that the MIT 
critical point in DMFT do have mean-fi eld critical behav- 
ior, confirming previous reports. ^ 9 * 11 ! From the Landau 
fit the critical temperature and coupling, T c and U c , was 
determined with high precision but with an accuracy lim- 
ited by finite size effects of the ED impurity solver. 

The study of the DMFT fixpoint surface using the 
phase space extension and the Jacobian has shown all the 
general features of the MIT. In principle it can be com- 
bined with numerical exact impurity solvers like Contin- 
uous Time Quantum MonteCarlo (CT-QMC) to remove 
convergence issues caused by fixpoint-bifurcations. This 
combination has the potential to further improve the ac- 
curacy in the location of the critical point. 

Our ED results for the critical point agrees with pre- 
vious ED calculations^! using the same number of bath 
sites, see Tab. |TJ Also the results of previous Hirsch-Fye 
QMC calculations^!! are compatible with our ED re- 
sults. The exact position of the critical point from IPT 
does not converge to the one of the other impurity solvers. 
This cannot to be expected since the expression for the 
self-energy is truncated at second-order in U. However 
our discretization proc edure yields results consistent with 
other IPT calculations P23 

Our results, showing that the fourth order Landau 
functional describes the critical properties of the DMFT 
solution when using ED, confirms the findings of Kotliar 
et. oZ.™ who showed, using Hirsch-Fye QMC and IPT, 
that the critical properties are not impurity solver de- 
pendent. 

Using the Landau functional we obtain the same uni- 
versal 1/3 critical exponents, in Eq. ( |36[ ), as initially re- 
ported for DMFT^ and later also found experimentally^! 
in Cr doped V2O3 in the U and pressure dependence 
respectively. We also present DMFT results on the 1/2 
critical exponent, Eq. ( 38 ) , along the first order transition 



line. This has been found experimentally in the tempera- 
ture dependence, D — D c ^ \T — T c \ 2 , of Cr doped V2O3, 
where the coupling to temperature in the first order term 
of the Landau functional vanishes, 70/ A) ~ 0. 

Regarding the minimal cusp singularity Landau model, 
Eq. ( 32 1 , it is noteworthy that the linear expansion of 



the a n (u, i) coefficients captures the critical behavior of 
the metallic and insulating solution of the DMFT data. 
However, an additional second order temperature term 

7o £ 2 must be added to correctly describe the unstable 
solution. Without this term, the Landau line U e (T) in 
Fig.[H]becomes a straight line, all isotherms of the Landau 
functional cross the same point (u — "fo/(3 t,d) — (0,0) 
(Fig. [9]) and the unstable central branch of the pitchfork 
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bifurcation in Fig. 10 becomes horizontal. 



V. CONCLUSION 

In this paper we have presented the phase space exten- 
sion algorithm for solving the single band DMFT fixpoint 
problem that is free from the numerical problems expe- 
rienced by other methods in the hysteresis region of the 
MIT. 

We have also performed an explicit calculation of the 
Jacobian of the DMFT fixpoint function and its corre- 
sponding eigenmodes. Using the critical eigenmode of 
the Jacobian we explained the critical slowing down and 
inability to find the thermodynamically unstable solu- 
tion using forward recursion. Moreover the instability of 
both forward recursion and Newton algorithms on the 
hysteresis boundaries has been explained in terms of one 
eigenmode going critical. 

The critical properties of the second order critical point 
of the MIT has been shown to be representable by a mean 
field Landau functional, giving the general properties of 
the DMFT fixpoint surface in (U, T, D) phase space in 
terms of a cusp singularity with a pitchfork bifurcation 



and saddle-node bifurcations on the hysteresis bound- 
aries. Experimentally the pressure driven MIT in Cr- 
doped V 2 O s show the same mean-field critical exponents 
for the second order critical point 

Finally we note that the phase space extension algo- 
rithm is general and can be combined with any impurity 
solver, and possibly also extended to study phase transi- 
tions in other systems such as the multi-band Hubbard 
model and the Hubbard-Holstein model. 
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